{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import rasterio\n",
    "#import matplotlib.pyplot as plt\n",
    "#%matplotlib inline\n",
    "from geopandas import GeoSeries\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_profile(ds,geom,dl,reverse=True):\n",
    "    \"\"\"\n",
    "    Uses an input rasterio dataset, a geometry in the form of Geoseries (polyline) and distance between points\n",
    "    to return the values of distance along the profile and the value from the dataset\n",
    "    \"\"\"\n",
    "    # reproject the shapefile to raster projection\n",
    "    geom = geom.to_crs(ds.crs)\n",
    "    #extracting point coordinates from the line at the specified distance\n",
    "    distance = 0\n",
    "    len = geom.iloc[0].length\n",
    "    while (distance<=(len)):\n",
    "        pt_geom = geom.interpolate(distance)\n",
    "        dicti = {'distance':distance,'geometry':pt_geom,'x':pt_geom.x.values,'y':pt_geom.y.values}\n",
    "        df = gpd.GeoDataFrame(dicti)\n",
    "        if distance==0:\n",
    "            df_f = df.copy()\n",
    "        else:\n",
    "            df_f = df_f.append(df)\n",
    "        distance = distance+dl\n",
    "    if reverse:\n",
    "        f_d = df_f['distance'].iloc[-1]\n",
    "        df_f['distance'] = np.abs(df_f['distance']-f_d)\n",
    "    #sampling the raster dataset values from the obtained point coordinates     \n",
    "    X = df_f.x.values\n",
    "    Y = df_f.y.values\n",
    "    xy = np.vstack((X,Y)).T\n",
    "    sampled_values = np.array(list(ds.sample(xy)))\n",
    "    return df_f['distance'].values, np.ma.fix_invalid(np.reshape(sampled_values,np.shape(sampled_values)[0]))\n",
    "    #return df_f['distance'].values, np.reshape(sampled_values,np.shape(sampled_values)[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_values(ds,geom):\n",
    "    \"\"\"\n",
    "    sample values from raster at the given ICESat-2 points\n",
    "    using nearest neighbour algorithm\n",
    "    \"\"\"\n",
    "    # reproject the shapefile to raster projection\n",
    "    geom = geom.to_crs(ds.crs)\n",
    "    #limit points within dataset bound\n",
    "    x_min,y_min,x_max,y_max = ds.bounds\n",
    "    #filter geom outside bounds \n",
    "    geom = geom.cx[x_min:x_max,y_min:y_max]\n",
    "    X = geom.geometry.x.values\n",
    "    Y = geom.geometry.y.values\n",
    "    xy = np.vstack((X,Y)).T\n",
    "    sampled_values = np.array(list(ds.sample(xy)))\n",
    "    return geom.x_atc.values, np.ma.fix_invalid(np.reshape(sampled_values,np.shape(sampled_values)[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
